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We study molecular para- hydrogen (P-H2) and ortho-deuterium (0-D2) in two dimensions and in 
the limit of zero temperature by means of the diffusion Monte Carlo method. We report energetic 
and structural properties of both systems like the total and kinetic energy per particle, radial pair 
distribution function, and Lindemann's ratio in the low pressure regime. By comparing the total 
energy per particle as a function of the density in liquid and solid P-H2 , we show that molecular 
para-hydrogen, and also ortho-deuterium, remain solid at zero temperature. Interestingly, we assess 
the quality of three different symmetrized trial wave functions, based on the Nosanow-Jastrow model, 
in the P-H2 solid film at the variational level. In particular, we analyze a new type of symmetrized 
trial wave function which has been used very recently to describe solid ^He and found that also 
characterizes hydrogen satisfactorily. With this wave function, we show that the one-body density 
matrix Qi{r) of solid P-H2 possesses off-diagonal long range order, with a condensate fraction that 
increases sizably in the negative pressure regime. 

PACS numbers: 61.50.Ah, 67.70.+n, 67.80.-s, 67.90.-fz 
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I. INTRODUCTION 



Quantum crystals like helium and hydrogen are in- 
triguing systems of fundamental physical interest. Due 
to the light mass of their constituents and relatively weak 
interparticle attraction, quantum solids exhibit large ki- 
netic energy and Lindemann's ratio even in the limit of 
zero temperature (T ~ 10° — 10"'^ K). In consequence, 
anharmonic effects and atomic quantum exchanges are 
of importance in this class of crystals. Moreover, in the 
last few years a series of ultra-low temperature experi- 
ments performed in solid ^He by different groups has led 
to a renewed interest on the possibility of superfluidity 
and/or Bose-Einstein condensation (BEC) in quantum 
solids lii^i^ii Essentially, these experiments analyze the 
quantum behavior of the helium crystal upon rotation or 
seek for some thermodynamic and/or structural anomaly 
signalizing a possible normal-to-superfluid phase transi- 
tion. Despite that most of the observations fairly agree 
in locating the onset of superfluidity (75 — 150 mK), there 
is a large dispersion in the value of the measured super- 
fluid fraction ps/p (~ 0.3 — 0.005 %). At present, there 
is lack of conclusive arguments for explaining these dis- 
crepancies but it is widely accepted that the purity of the 
sample and the presence of crystalline defects play on it 
a relevant rolei^ On the other hand, there is overall 
agreement among microscopic full quantum calculations 
in practically ruling out superfluidity in the perfect (free 
of defects) bulk configuration. It is worth noticing that 
usual techniques devised to study classical crystals (that 
is crystals composed of heavier elements and with larger 
cohesive energies), like for instance harmonic based ap- 
proaches, are not longer suitable for quantum solids and 



calculations on them are in most cases challengingiii^ 

In the present work, we present a theoretical study 
of two-dimensional (2D) molecular para-hydrogen (p- 
H2) and ortho-deuterium (0-D2) at zero temperature by 
means of the diffusion Monte Carlo (DMC) method^ d^d^ 
and the semi-empirical radial pair interaction due to Sil- 
vera and Goldman. ^^ Hydrogen is a very interesting and 
challenging system which has been investigated very in- 
tensively during the last half century. As a matter of 
fact, hydrogen is the most abundant element in the uni- 
verse and from a technological point of view it is con- 
sidered among the most promising green combustibles of 
the near future. Very interestingly, hydrogen has been 
predicted to exhibit a new state of matter at very high 
pressures {P ~ 400 GPa) in which superfluidity and su- 
perconductivity might coexist ji^ii^ 

In this work, we restrict our analysis of molecular hy- 
drogen and deuterium to the low pressure regime (P ~ 0) 
and zero temperature. Contrarily to what occurs in he- 
lium, molecular hydrogen freezes at a temperature of 
13.96 K in spite of its lighter mass, given that the interac- 
tions between particles are more attractive (the minimum 
of the interaction between hydrogen molecules amounts 
to ^ —37 K while in helium is ~ — 10 K). One of our moti- 
vations for carrying out the present study was to unravel 
whether liquid P-H2 could be stabilized or not at zero 
temperature by reducing the dimensionality with respect 
to the bulk. This possibility appears to be very appeal- 
ing since it would provide a chance for superfluidity and 
Bose-Einstein condensation (BEC) to be observed in a 
quantum liquid different from helium. In fact, P-II2 in 
one dimension and inside a carbon nanotube has already 
been studied in the zero-temperature limit and predicted 
to be liquid at its equilibrium densityJ^ Also, small drops 



with a number of molecules A'^ < 26 present superfluid 
character J^^y^'— On the other hand, it has been reported 
recently an experiment performed on molecular ortho- 
deuterium pre-plated on krypton at very low tempera- 
tures (T ~ 1 K) in which it is claimed evidence for the 
existence of a reentrant 0-D2 liquid phasej^^ As it will be 
presented in short, our results show that no first order 
solid-liquid phase transition occurs in two-dimensional 
H2 at zero temperature. On account of this result, we 
straightforwardly reject this possibility also for 0-D2 since 
deuterium molecules are heavier and their intermolecu- 
lar interactions are considered equal to the H2-H2 ones. 
Consequently, most of the effort done in this work has 
been devoted to achieve an accurate description of the 
ground state of two-dimensional solid P-H2 and 0-D2 . 

In this work, we calculate by means of the diffu- 
sion Monte Carlo method (DMC) some energetic and 
structural properties of both hydrogen and deuterium 
films near equilibrium. Quantities like the kinetic en- 
ergy per particle and Lindemann's ratio have been com- 
puted within the pure estimator approacb^^^i^ in or- 
der to remove any possible bias coming out from the 
trial wave function used for importance sampling. In 
this way, we quote quantum isotopic effects in hydrogen 
directly and only within the statistical uncertainty. Pre- 
vious to the DMC results, we present a variational Monte 
Carlo (VMC) study of the P-II2 crystal in which we have 
tested the quality of several symmetrized and unsym- 
metrized trial wave functions. With this analysis, we de- 
termine the effect of symmetrization on the total energy 
and the relevance of molecule exchanges along the sim- 
ulation. Moreover, we analyze which symmetrized wave 
functions can be implemented in DMC to the end of esti- 
mating the possible superfluidity and Bose-Einstein con- 
densation (BEC) of the solid with simultaneous accurate 
description of their energetic and structural properties. 
In particular, we have studied in detail a symmetrized 
trial wave function, named ipjQ in this work (the nota- 
tion will become clear later), which has been proposed 
and used recently to study bulk solid '^Hc.-^ Here, we 
find that iJjjq also characterizes solid hydrogen in two 
dimensions accurately. 

Interestingly, we assess the behavior of the one-body 
density matrix gi{r) of P-II2 with density by means of the 
symmetric trial wave function V'/g- ^^ ^^^ cases, a very 
small condensate fraction uq is observed. For densities 
below the equilibrium one, and near the spinodal point, 
a significant increase of no is observed pointing to the 
emergence of a finite superfiuid density. 

This paper is organized as follows. In Sec. |TT1 we 
present a brief description of the semiempirical pair po- 
tential and techniques used throughout this work. Next, 
in Sec. lIIII and llVl we report our variational and diffusion 
Monte Carlo results for P-H2 and 0-D2 in two dimensions, 
respectively. Sec. |V]is devoted to the examination of the 
one-body density matrix gi (r) obtained with 4'jg ^^^ ^^^ 
dependence with the density. Finally, in the Sec. I VII we 
summarize the main results presented in this work. 



II. MOLECULAR INTERACTION AND 
METHOD 

The H2 (D2) molecule is composed of two hydrogen 
(deuterium) atoms linked by a covalent bond, which 
in the para-hydrogen (ortho-deuterium) state possesses 
spherical symmetry (total angular momentum zero). The 
energy scale involved in electronic excitations (~ 10^ K) 
is orders of magnitude larger than the intermolecular one 
(~ 10^ K), thus to model the H2 - H2 (or D2 - D2) in- 
teraction by means of a radial pair-potential and consider 
the molecules as point-like turns out to be justified upon 
the condition of low or moderate pressures. In this work, 
we have adopted the well-known and commonly used 
semiempirical Silvera-Goldman pair potential/^ This po- 
tential has proved to perform excellently at low temper- 
ature and at the pressure regimes in which we are inter- 
ested. 

The ground state of para-H2 and ortho-D2 is de- 
termined using the DMC method. DMC is a 
zero-temperature method which provides the exact 
ground-state energy of many-boson interacting systems 
within some statistical uncertaintyi^ i^°'^^ This tech- 
nique is based on a short-time approximation for the 
Green's function corresponding to the imaginary time(T)- 
dependent Schrodinger equation. Despite this method 
is algorithmically simpler than domain Green's function 
Monte Carlofiii^ it presents some (Ar)" bias coming 
from the factorization of the imaginary time propaga- 
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Our implementation of DMC is quadratic 
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hence the control of the time-step bias is efficiently con- 
trolled given that the required At — > extrapolation is 
nearly eliminated by choosing a sufficiently small time 
step. The Hamiltonian H, describing our system is 
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and the corresponding Schrodinger equation in imaginary 
time {it = r). 



_c)^(R,r) 



= (H-i?)*(R,r), 



(2) 



with E an arbitrary constant. Equation ([2|) can be 
formally solved by expanding the solution \['(R, t) in 
the basis set of the energy eigenfunctions {^n} (R- = 
{ri, r2, . . . , rjv}). It turns out that ^^(R, r) tends to the 
ground-state wave function $0 of the system for an infi- 
nite imaginary time as well as the expected value of the 
Hamiltonian tends to the ground-state value Eq. The 
hermiticity of the Hamiltonian guarantees the equality 
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where ipx is a convenient trial wave function. As a con- 
sequence, the ground-state energy of the system can be 



computed by calculating the integral 



{ii)DMC - lim / El (R) / (R, r) dR , 
■^^°° Jv 



(4) 



where / (R, r) = * (R, r) Vt (R), and El (R) is the lo- 
cal energy defined as £'l(R) = Uipr (R) /V^t (R)- The 
introduction of V't (R) in / (R, r) is known as impor- 
tance sampling and its use is important to reduce the 
variance of Eq. (j4|) to a manageable level (for instance, 
by imposing f/r (R) ~ when rij is smaller than the core 
of the pair interaction). 

In this work, all the operators diagonal in real-space 
which do not commute with the Hamiltonian, that is 
[H, O] ^ 0, have been sampled using the pure estimator 
technique based on forward walking i22i2iiS^ Essentially, 
with this method the possible bias introduced by '0T in 
the mixed estimator ('&o|O|'0t) are removed by proper 
weighting of the configurations generated along the sim- 
ulation. 



III. MOLECULAR PARA-HYDROGEN 

A. Variational Monte Carlo results 

In this section, we present a variational Monte Carlo 
(VMC) study of two-dimensional P-H2 which provides 
us with the most convenient trial wave function (twf) 
to be used in subsequent DMC calculations and valu- 
able physical insight on the system itself. In brief, the 
VMC method relies on the variational principle which 
states that given a Hamiltonian the energy difference 
El — Eq averaged over the probability density distri- 
bution IV'Tp is always positive and it decreases as the 
overlapping between ipx and the true ground-state wave 
function increases {El and E^ are the local energy de- 
fined in the previous section and the ground-state energy, 
respectively). In the present work, the main variational 
effort has been devoted to achieve an accurate descrip- 
tion of the 2D solid phase. We have checked by means of 
VMC and DMC that in both P-H2 and 0-D2 systems the 
triangular configuration is the stable one at all the stud- 
ied densities. The energies reported in this section have 
been calculated at the density p = 0.060 A^^ , which 
corresponds to the variational equilibrium density of the 
solid. 

In order to determine the nature of the ground-state of 
the system we have also carried out simulations for the 
liquid phase. In this case, the trial wave function is of 
Jastrow type, 
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■0j(ri,r2,...,rjv) ^H^^C^ 
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where the two-body factors f2 account for the molecular 
correlations arising in the system due to pair interactions. 
These two-body correlation factors have been chosen of 



McMillan form, f2 = e~2(,7J ^ and as best value of the 
variational parameter b we obtain 3.70 A . 

For the solid phase, an additional one-body factor is 
introduced ad hoc in the trial wave function to the end of 
reproducing the periodic order of the system and so mak- 
ing the sampling over the space of configurations more ef- 
ficient (Nosanow- Jastrow model). Such one-body factor 
consists in a productory of localizing functions centered 
on the positions that define the perfect crystal configu- 
ration (sites), given by the family of vectors {Ri} 

N 

ipNj{ri,r2,...,rN) = i'jYigiUrt -Ri|) ■ (6) 

i 

We have explored two different trial wave functions based 
on ipNj , each one consisting in a different choice for 
gi (r) . The first model corresponds to the standardly used 
Gaussian function, 



gcir) =exp(^-^r2j 



while the second is a Fade function, defined as 

apr^ 



.(r) = exp - 
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(8) 



where ac ap and cp are variational parameters to be 
optimized. Even though gci''') and gp{r) are analyti- 
cally quite similar, the asymptote of go can be chose 
so as to decay to zero less abruptly than that of gair) 
(see Fig. [IJ. This feature can be used in the simulations 
for increasing somewhat the degree of delocalization of 
the molecules at the expense, however, of an increase of 
the kinetic energy within the surroundings of the equi- 
librium positions (where gp{r) varies more rapidly than 
gG(?') )■ In Table I , we report the best energies obtained 
in the optimization process of ipNj with Gaussian and 
Fade functions for gi(r) as well as the optimal set of 
parameters. As one can see, in both cases the lowest 
energy obtained is —21.3(1) K . It is worthwhile notic- 
ing that when the asymptotes of the Fade factors are 
widened (that is, the value of cp is increased), or equiv- 
alently, when the molecules are left to move more freely 
around the equilibrium positions, the energy of the sys- 
tem increases. Given the variational equivalency between 
gG(r) and gp(r) , we opt for V'a'J with Gaussian factors 
and optimal parameters ba — 3.45 Aand gq — 0.67A~^ 
in our subsequent DMC calculations. 

It is well known that iJjnj is not symmetric under the 
permutation of particles (that is, '0Arj(ri, r2, . . . ,r7v) j^ 
V'w j(i"2, Ti, . . . , rjv)). This property is manifestly incor- 
rect for a system of indistinguishable bosons. Neverthe- 
less, the use of the Nosanow- Jastrow model is widely 
spread within the field of microscopic calculations since 
it is assumed that the effect of symmetrization on the 
total and partial energies of quantum solids is practically 
negligible. In fact, we will show in brief that this is also 
the case for two-dimensional hydrogen. To this end, we 




FIG. 1; Optimized Gaussian and Fade functions (solid and 
dashed line, respectively) for solid P-H2 at the density p — 
0.060A~'^ . The decay of gp (r) to zero is smoother than that of 
gcir), and in the region around the origin gp{r) is narrower. 



have tested the quality of several symmetrized trial wave 
functions in the variational study of 2D solid P-H2. As 
mentioned in the Introduction, our motivation for this 
analysis is twofold: estimate the influence in the energy 
of molecular quantum exchanges, and chiefly, study the 
possibility of its use as importance sampling in DMC to 
determine the possible existence of Bosc-Einstein con- 
densation (BEC) and superfluidity. 

We have studied three different symmetrized trial wave 
functions. The first model consists in a permanent of 
monoparticular functions containing the A^! possible per- 
mutations {P} of the N particles among the different 
lattice sites, expressed as 



N 



V'ji(ri,r2,. 



> rw) = V'J X! n gi ('■■' ^ ^P' 

{P}^=1 



(9) 



Due to the algebraic difflculties arising in the implemen- 
tation of permanents (contrarily to what occurs with de- 
terminants), the sampling of ipjj^ must be divided into 
two different parts, one performed in the space of spa- 
tial configurations and the other in the space of per- 
mutationsi^^ The acceptance probability for a proposed 
change of position of the particle labeled i, r^ ^ rj-, cor- 
responds to 

.^minfl,^^(-^g^^<-^-)gf-^-)V (10) 
V V',/(r)2gi(r, - ROgi(r. - Rp.) J' ^ ' 

where the subindex Pi can take any of the A^ possible 
lattice sites. On the other side, the acceptance proba- 
bility for a proposed site permutation between the i and 
the j particles, Rpi ■s-^ Rpji is 



Q2 



^ gi(rj -RpOgi(r, -Rpj) 
'gi(rj -Rpj)gi(rj -Rp.O 



(11) 



Notice that permutations involving more than two parti- 
cles are not sampled since the acceptance level for swap 
permutations is already extremely low. 



The optimal results obtained with ipjj^, using Gaus- 
sian and Pade gi(r) functions, are reported in Table I . 
By comparing the variational energies obtained with V'Af J 
and ipjj^ , we show that symmetrizing tpNj with the above 
prescription has not appreciable effects on the total en- 
ergy of P-H2. Nevertheless, it must be said that one 
should not draw other conclusive statements about the 
effects of a full symmetrization just based on the approx- 
imation of the permanent by a reduced sampling in the 
permutation space of the type pT|) . In fact, the ac- 
ceptance rate of permutations is so low (column Q2 in 
Table I ) that sampling ipfj^ efficiently turns out to be 
quite challenging. 

The second model of symmetrized trial wave function 
tpj'i^ consists in a productory of sums in the form 



Af 
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i^jl (i-i, r2, . . . , Tat) = V'j n Zl Si (r^ - Rj 

(12) 

This trial wave function has been proposed very recently 
by Zhai and Wu'^^ and has been suggested to be of possi- 
ble relevance for the study of the supersolid. In fact, tpjl 
avoids any explicit sampling in permutation space hence 
turns out to be well-suited for being used as importance 
sampling in DMC simulations. However, as one can see 
in Table I the best variational energy obtained with this 
model is sizably larger than the ones obtained with V'atj 
and V'fi; in both Gaussian and Fade cases. In fact, the 
variational energy obtained with -^j^ is very similar to 
the one calculated with tpj for the liquid (~ —17.4 K). 
In doing the simulation with this wave function it is ob- 
served that particles diffuse excessively within the con- 
tainer giving place to glassy-like configurations; we have 
checked this feature by monitoring the radial pair distri- 
bution function (see Fig. [2]) and mean squared displace- 
ment (which grows steadily with time). The reason for 
this excessive atomic diffusion is that the way in which 
the one-body factor is symmetrized in tjjj^ does not pe- 
nalize multiple occupation of a same lattice site. This 
feature will be illustrated in short by means of a simple 
example involving two particles moving in one dimen- 
sion. Moreover, if the width of gi(r) is narrowed in order 
to avoid such unrealistic molecular diffusion, the total 
energy of the system is worsened because of the rapid 
increase of kinetic energy. 

The third type of symmetrized trial wave function 
reads 
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V'jG(r 1,1-2, 



,T^N 



^jU 



;gi(r,-R,) , (13) 
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and it is also straightforward to implement in DMC 
codes. This type of trial wave function has been pro- 
posed very recently by Cazorla et alM' and has been used 
to set an upper bound of 10~^ for the superfluid fraction 
of perfect crystalline bulk ''He at zero temperature. tpjQ 
and ipjQ look similar, the difference being that the pro- 
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TABLE I: Optimal variational total energy per particle obtained at the density p= O.O6OA ^ with different trial wave func- 
tion models. Values appearing on the left(right)-side of the table correspond to one-body factors gi(r) adopted in the form 

&G{r) (gp(r)) . 
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FIG. 2: Variational radial pair distribution function g{r) 
of two-dimensional molecular hydrogen at the density p — 
0.060A"^ obtained with twf ip% and V'fS • 
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FIG. 3: Squared Vjc and t/^jQ (with t/jj = 1) in the sim- 
ple case of two particles moving in one dimension and sites 
separated by one arbitrary unity. 



ductory and summatory in tp^Q run over sites and parti- 
cles, respectively, while in ipjQ is the other way around. 
However, ipjQ and ipjQ lead to completely different vari- 
ational energies (see Table I ). The best variational re- 
sult obtained with ^jq amounts to —20.4(1) K, which 
is only 0.9 K larger than the one calculated with V'jvj 
or ipjj^ but 2.5 K smaller than the one corresponding to 
tpjQ- In this case, the radial pair distribution function 
and mean squared displacement follow typical solid-like 
patterns (see Figl^]). 

Contrarily to what occurs with tpjQ, the multiple oc- 
cupation of a same site is now penalized by the wave 



function and hence crystal order is sustained. To the end 
of illustrating this feature, which appears to be the main 
difference between tpjQ and ipjC' ^^ have analyzed the 
simple case of two particles moving in a one-dimensional 
lattice. For the sake of simplicity, we have assumed that 
the distance between the equilibrium positions of the 
particles is one, that the parameter entering the Gaus- 
sian factors in Eq.(I2]) and P^ is a = 1/2, in arbitrary 
units (a.u.), and that the Jastrow factor is switched off 
{tpj = 1). The value of the squared wave function for tpjQ 
and ipjQ, \ipsoi\'^, obtained by keeping fixed one of the 
particles in the site located at the origin and then moving 
the other particle towards it, is plotted in Fig. [3] for the 
interval < a; < 1 . As one observes in there, the value of 
4>jQ at X = 1 and (that is, each particle is placed over 
one site or both are at the same position, respectively) 
is the same, whereas i^jQix = 1) > ipjai^ — 0)- This 
effect is what we have previously refereed to as "penal- 
ized by the trial wave function" . It is also noted that the 
value of tpjG ^^ niaximum at half the way between and 
1, not so for ifijQ, hence ipj^ will always promote larger 
diffusion of the molecules. 



B. Diffusion Monte Carlo results 

We have studied the energetic and structural prop- 
erties of P-H2 using the DMC method and ipNj & as 
trial wave function. We have verified that the DMC en- 
ergy and diagonal properties obtained with ipNj are sta- 
tistically indistinguishable from the ones obtained using 
the symmetric wave function tpjG P^p . The results pre- 
sented in this section have been obtained for a rectangu- 
lar plane containing 90 particles with periodic boundary 
conditions in the two spatial directions. Internal param- 
eters of the simulations, namely the averaged population 
of walkers and time step, are 250 and 5 ■ 10"'* K"-'^, re- 
spectively; these parameters have been adjusted in order 
to reduce any possible bias to the level of the statistical 
uncertainty (~ 0.05 K). 

In Tabic II , we report the total ground-state energy 
per particle, E/N, corresponding to 2D solid P-H2 at 
some densities. The pure (unbiased) estimation of the po- 
tential, V/N, and kinetic energies, T/N, are also quoted 
therein. The energy results have been corrected for the 
finite size of the simulation plane by assuming the radial 



p(A-^) E/N (K) V/N (K) T/N (K)" 

0.053 -19.42(2) -35.73(3) 16.30(3) 

0.060 -22.21(2) -43.67(3) 21.46(3) 

0.065 -23.27(2) -49.23(4) 25.96(4) 

0.067 -23.42(2) -51.68(4) 28.26(4) 

0.070 -23.19(2) -54.83(4) 31.64(4) 

0.076 -21.30(2) -59.22(5) 37.92(5) 

0.083 -14.23(2) -62.40(7) 48.17(7) 



TABLE II: Ground-state energy E/N, potential energy V/N, 
and kinetic energy T/N, per particle of solid 2D P-H2. Poten- 
tial and kinetic energies are obtained with pure estimators. 



pair distribution function g{r) to be one beyond the dis- 
tance i?max = L/2, with L being the size of the plane. 
Assuming g{r) ~ 1 beyond i?niax could seem a crude 
approximation for crystals since this function shows pe- 
riodic structure (see for instance Fig. [5|). However, the 
periodic oscillations of g{r) around unity might suggest 
that in average this approximation is essentially correct. 
In order to test the reliability of this finite size correc- 
tion we have calculated the total energy per particle in 
a plane containing 90 , 120 and 168 molecules at the 
density p = 0.0597 A^^ ; we obtain E/N = -22.19(2), 
—22.16(2) and —22.15(2) K, respectively, thus achieve- 
ment of convergence within the present statistical uncer- 
tainty is proved. 

The energy per particle corresponding to liquid and 
solid 2D P-H2 at zero temperature is plotted in Fig.|4]as a 
function of the density. The simulation of the metastable 
liquid phase uses a Jastrow wave function ■0,/ ^ as im- 
portance sampling. The lines in Fig. |3] correspond to 
polynomial fits to our results in the form 



E/N = e{p)^eo + B 
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The pressure, compressibility and speed of sound (aver- 
aged for all the directions) are then easily derived from 
Ea.([T4)) through the expressions. 
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P \dP 



f 1 



\mKp 



(15) 



(16) 



(17) 



The optimal value of the parameters for the solid phase 
are eg = -23.453(3) K, po = 0.0673(2) A-^, B = 
121(2) K, and C = 152(8) K, where eo and po are 
the equilibrium energy per particle and density, respec- 
tively. According to these figures, the compressibil- 
ity and speed of sound at the equilibrium density are 



0.055 0.06 0.065 0.07 0.075 

P(A-=) 

FIG. 4: Total ground-state energy per particle of liquid (dot- 
ted line) and solid (dashed line) 2D P-H2 at zero temperature. 
The lines correspond to polynomial fits of our results (empty 
and filled circles); the statistical errors bars are smaller than 
the symbol size. 



k(Po) = 0.0615(8) AVK and c(po) = 998.6(1) m/s, the 
numbers quoted within parentheses being the statistical 
errors. The equation of state of the hquid phase is also 
well described by the polynomial form (fTi|) with opti- 
mal parameters eo = -21.43(2) K, po = 0.0633(3) A'^, 
B = 75(7) K, and C = 69(9) K . 

Another magnitude of interest in the study of bulk 
systems is the spinodal density ps ■ ps sets the limit 
for the system to remain in a homogeneous phase since 
at this density the compressibility grows to infinite (or 
equivalently, the speed of sound becomes zero); in case 
of going below this point (p < ps) the system breaks 
down into clusters. According to our DMC calculations, 
this low-limit density amounts to ps = 0.0548(1)A^^ in 
solid P-H2. 

A glance at Fig. 2] shows that the solid phase is the 
stable one overall the regime of positive pressures. Nev- 
ertheless, by looking at our results one could suggest a 
first order liquid-solid phase transition occurring at nega- 
tive pressures, where the two equations of state cross each 
other. Needless to be said, that this possibility deserves 
detailed exploration since it could provide a chance for 
superfluidity to be observed in a quantum liquid different 
from helium. Aimed at this, we have simulated 2D liquid 
P-H2 down to densities of 0.039A~^ (the spinodal den- 
sity of the liquid is ps = 0.0519(1)A^^) and searched for 
that transition by means of the Maxwell double-tangent 
construction. Our results show that a liquid-solid transi- 
tion is not possible within the range set by the spinodal 
densities, and thus the possibility of liquid P-H2 in two 
dimensions must be rejected. 

Our results for the equation of state of P-H2 can be 
compared with two previous PIMC studies carried out 
on the same system. In Ref 28, Gordillo and Ceperley 
obtained p = 0.064 A^^ for the equilibrium density of 
2D solid P-II2 at T = 1 K; the authors of that work 
reported a figure with the energy per particle as a func- 
tion of the density, and the minimum of the curve is 
located at ~ —22.0 K. A more systematic analysis of the 



2.5 




2 


1 


1 


hh'h-fv 


0.5 


1 |/ \W v ^■■■'' 



r(A) 



FIG. 5: Radial pair distribution functions of 2D solid 0-D2 at 
the equilibrium density po — 0.078 A~^ (dashed line), and of 
P-H2 at the density 0.068 A~^ (solid line). 



P(A-^) 


7H2 


C(10) 


C(01) 


(u^)(A^) 


0.058 


0.212(1) 


0.00(2) 


0.00(2) 


0.19(1) 


0.060 


0.197(1) 


0.02(2) 


0.00(2) 


0.15(1) 


0.065 


0.183(1) 


0.02(2) 


0.02(3) 


0.12(1) 


0.067 


0.178(1) 


-0.01(1) 


-0.02(2) 


0.11(1) 


0.070 


0.170(1) 


-0.02(1) 


0.00(2) 


0.09(1) 


0.076 


0.158(1) 


0.00(2) 


0.01(1) 


0.07(1) 


0.083 


0.146(1) 


0.01(2) 


0.01(1) 


0.06(1) 



TABLE III: Lindemann's ratio 7H2, kurtosis C, and mean 
squared displacement (u'^) of 2D solid P-H2 at different den- 
sities close to equilibrium (pure estimations). 



same system was performed later on by Boninsegnij^ 
In that work, the total energy per particle and chemical 
potential are calculated at several densities and within 
the temperature range T = 1 — 8 K. Subsequently, an 
extrapolation of the low temperature results to absolute 
zero was performed, leading to Pq^^'^ = 0.0668(5) A~^, 
gPiMC ^ _23.25(5) K and pf^^ = 0.0585(10) A~^. We 
note that the agreement between those zero-temperature 
extrapolated PIMC values and our DMC results is fairly 
good, specially in the case of the equilibrium density po- 
We have analyzed the structure of the 2D solid by cal- 
culating the radial pair distribution function g(r). 



sW = 



iV - 1 / I* (ri, ri + r, . . . , Tat) l^dndrg 



. dr 



N 



p /|*(ri,r2,...,rAr) |2dri...(ir 

and the Lindemann's ratio 7H2 , 



N 



7 



\ 



1 ^ 
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(19) 



where a is the distance between nearest neighbors in the 
perfect crystalline configuration. In Fig. [51 we plot g{r) 
at the density 0.068A~^ which, as it is expected in crys- 
tals, exhibits a pattern of periodic order. At low temper- 
atures, the Lindemann's ratio 7 around the equilibrium 
density tends to zero in classical solids while in quantum 



crystals it is finite due to the zero-point motion of parti- 
cles, hence this quantity is regarded as a good quantum 
indicator. Furthermore, the Lindemann's ratio (or equiv- 
alently, the mean squared displacement) is related to the 
Debye- Waller factor Mq , which describes the attenua- 
tion of the emergent radiation in coherent scattering ex- 
periments according to the formula I{Q,T) oc e'"^^"'^''^ 
(where I {Q,T) is the intensity of the outgoing radiation 
scattered by the target and Q is the modulus of the trans- 
fer wave vector). By means of a cumulant expansion, the 
Debye- Waller factor can be expressed as 



2M, 



Q 



,)Q' 



12 



{{u%)^Hul)^)Q^ + 0{Q'), (20) 



where {u%) is the mean squared displacement along the 

direction Q. It is easy to see that when the distribu- 
tion of particles around the equilibrium positions is well- 
described by a Gaussian function, the quantity within 
parentheses in the second right-term of Eq. (|^D)) , known 
as kurtosis Qq , vanishes. In such a case, the Debye- 
Waller factor reduces to the simple formula 2Mq = 
{uq)Q^ . In Table III , we report the Lindemann's 
ratio, kurtosis and mean squared displacement of two- 
dimensional P-II2 at different densities. We have calcu- 
lated (uq) and Cq along two orthogonal directions and 
not found appreciable differences in the results. More- 
over, the kurtosis is null in all the studied cases. Conse- 
quently, we may conclude that the distribution of hy- 
drogen molecules around the equilibrium positions is 
isotropic and can be accurately reproduced by a Gaus- 
sian, contrarily to what it is found to occur in ■*He^ 
Regarding the value of 7, it can be said that solid H2 is 
less quantum than "^He since 7H2 ~ 0.18 at po whereas in 
solid helium 7hc ^ 0.24 near melting.^ Also it is worth 
noticing that the trend of 7H2 is to increase with de- 
creasing density, therefore quantum exchange effects in 
the crystal would become of greater relevance at small 
densities. 



IV. MOLECULAR ORTHO-DEUTERIUM 

The ground-state properties of 0-D2 (with total an- 
gular momentum zero) have been also studied using the 
DMC method and the same radial pair potential (Silvera- 
Goldman) than in P-H2. The larger mass of D2 makes 
one to expect that two-dimensional bulk D2 is solid at 
zero temperature, so in this case we have restricted our 
study to the solid phase. 

In our simulations, the equilibrium positions of the o- 
D2 molecules are arranged according to a triangular lat- 
tice and the particles are assumed point-like. In this case, 
we use the trial wave function 



■^Arj(i'i:i"2,...,rjv) = n^ ' 



N N 

|(fc/n,) 



n 



e 2 



(|r,-R, 



(21) 
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FIG. 6: Ground-state energy per particle of 2D solid 0-D2 
(solid line and filled circles), and 2D solid P-H2 (dotted line 
and empty triangles) which is shown for comparison. 



P (A-) 


E/N (K) 


V/N (K) 


T/N (K) 


0.046 


-22.40(2) 


-29.84(3) 


7.44(3) 


0.053 


-27.97(1) 


-37.70(3) 


9.73(3) 


0.060 


-33.57(1) 


-46.38(4) 


12.81(4) 


0.069 


-39.41(1) 


-57.52(3) 


18.11(3) 


0.076 


-42.12(1) 


-66.46(6) 


24.34(6) 


0.084 


-40.80(1) 


-72.59(6) 


31.79(6) 


0.094 


-29.32(2) 


-73.18(8) 


43.86(8) 



TABLE IV: Ground-state total and pure potential and ki- 
netic energies per particle of 2D 0-D2 at several densities. 



which differs slightly from ipNj in Eq. ([6|) ( now the pair 
correlation factors £2 depend on the extra variational pa- 
rameter c ) . The variational parameters in Eq. (|2ip have 
been optimized using VMC; the best values are 6 = 3.32 
A , c = 7, and a = 0.67 A'^ . All the DMC simulations 
have been performed in a rectangular plane containing 
120 particles and applying periodic boundary conditions. 
The target averaged population of walkers, n^, is 250 
and the time step, Ar, 5 • lO^"' K"-'^ . Finite size effects 
have been corrected with the same approach than used 
for hydrogen (see Sec. lIIIB")) . 

In Fig. [5] , we plot the total ground-state energy per 
0-D2 molecule as a function of the density; the solid line 
represents the best fit to our data following the poly- 
nomial function expressed in Eq. P^ . The best value 
of the parameters are B = 241(3) K, C = 324(10) K, 
eo = -42.305(5) K and po = 0.0785(2) A^^ _ ^Ya^^v lead 
to a spinodal density ps — 0.0641(2) A^^ . By com- 
paring with P-H2 , we show that 0-D2 is denser at equi- 
librium and appreciably more bounded (the total energy 
decreases substantially). The heavier mass of the 0-D2 
molecules, makes the solid to reduce its kinetic energy 
and mean squared displacement at any density (see Ta- 
bles II , III , IV and V ), thus allowing the system to 
increase its equilibrium density in order to take advan- 
tage of the attractive interparticle interaction. 

Concerning the structural properties of 2D solid 0-D2 , 
we have calculated the radial pair distribution function 



P(A-) 


7D2 


C(10) 


C{01) 


0.053 


0.204(1) 


0.00(1) 


-0.06(2) 


0.060 


0.187(1) 


0.01(1) 


0.00(1) 


0.069 


0.160(1) 


0.00(2) 


-0.01(1) 


0.073 


0.149(1) 


-0.01(1) 


0.01(1) 


0.078 


0.139(1) 


0.00(1) 


0.01(1) 


0.080 


0.135(1) 


-0.02(1) 


0.00(1) 


0.088 


0.124(1) 


0.00(1) 


-0.03(1) 


0.094 


0.117(1) 


-0.02(1) 


0.02(1) 



TABLE V: Lindemann's ratio 7D2! ^-nd kurtosis C, of two- 
dimensional 0-D2 at different densities near the equilibrium. 



at the equilibrium density (see Fig. [5] ), and the Lin- 
demann's ratio and kurtosis at different points (see Ta- 
ble V ). As it is shown in Fig. [5] , the peaks of the radial 
pair distribution function of 2D 0-D2 are sharper and 
somewhat closer than in molecular hydrogen at equilib- 
rium since the density and degree of localization of the 
particles are larger in the first case. This statement is 
also corroborated by the results contained in Table V , 
where the Lindemann's ratio is invariably some tenths 
smaller than in P-H2 at the same density. As a matter 
of comparison, the Lindemann's ratio of 0-D2 at equilib- 
rium is about 1.3 times smaller than that of P-H2 and 1.7 
than in two-dimensional "^He . Furthermore, as it could 
be expected from our previous study of 2D P-H2 (Table 
IIIIBp , the kurtosis in two-dimensional 0-D2 is practically 
null in both the two orthogonal directions for which it has 
been calculated. 



V. ONE-BODY DENSITY MATRIX AND 
OFF-DIAGONAL LONG RANGE ORDER 

A fundamental function in the study of quantum sys- 
tems is the one-body density matrix pi(r,r'), defined as 



gi(r,r') = ($o|^^(r)^(r')|*o) , 



(22) 



where V'(r') and ?A^(r) are, respectively, the field opera- 
tors which destroy a particle from position r' and create 
another at position r, and $0 is the ground-state wave 
function. In particular, a finite value for \iuYr—>ooQi{f) 
proves the existence of off-diagonal long range order 
(ODLRO) in the system, the measure being the conden- 
sate fraction no- In quantum Monte Carlo, the one-body 
density matrix can be estimated by averaging the coordi- 
nate operator Vrlri + r, r2, . . . ,rAr)/V'T(ri, r2, . . . ,rAr) . 
Here, we use extrapolated estimators for qi since the pure 
estimation relying on forward walking is only applicable 
to diagonal operators. Moreover, in order to get consis- 
tent results we have required that the two extrapolated 
estimators of the same accuracy, i.e. , Qi{r) = 2g^'^^{r) — 

gi"'^{r) and gi{r) — (£'i"'^(f)) /gi"'^{r) (where mix 
means obtained with DMC and var with VMC), coin- 
cide within the present statistical uncertainty. 
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FIG. 7: One-body density matrix of P-H2 Qi{r) obtained 
using as importance sampling ipNj and ipjc ^^ the density 
p — 0.060 A~^ . The sohd line corresponds to the Gaussian 
function which best fits to the Nosanow-Jastrow result. 
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FIG. 8: Top: One-body density matrix of P-H2 at density 
p = 0.060 A~ obtained using importance sampling with 
the Nosanow-Jastrow and ipjG trial wave functions. Bottom: 
One-body density matrix of P-H2 obtained with the tpjo trial 
wave function at several densities. Densities are in units of 

A-2. 



In Fig. [71 we compare DMC results of the one-body 
density matrix at the density p = 0.060 A~^ and within 
the distance range < r < 5.0 A , obtained with both 
4'Nj (unsymmetrized) and tpjo (symmetrized) trial wave 
functions. As one can observe therein, the series of points 
obtained with both twfs are compatible in the full de- 
picted range. In the same graph, we also enclose the 
Gaussian curve G{r) = e"*"" which best fits to the re- 
sult obtained with ipNj ; we find in this case that the 



optimal value of the parameter b is 0.400(6) A~^ . In 
order to test the quality of this fit (which in the reduced 
chi-squared test leads to the value 1.59), we have calcu- 
lated the atomic kinetic energy of two-dimensional P-II2 
through the formula 



T/N 



2mH 



-V^giir) 



(23) 



r=0 



but assuming G{r) instead of gi{r). In fact, it may be 
shown that Eq. (j23p derives from the second moment of 
the momentum distribution n(fc). 



T/N 



1 



2toh2 (27r)" 



dk k^ n{k) . 



(24) 



Proceeding in this way, we obtain T/N = 19.40(30) K 
which does not agree satisfactorily with the correspond- 
ing pure (mixed) estimation 21.46(3) (20.75(2)) K . Very 
interestingly. Withers and Glyde have recently shown by 
means of simple models that anharmonic and/or particle- 
exchange effects in quantum solids may cause the mo- 
mentum distribution n(k), or equivalently gi{r), to de- 
viate significantly from a Gaussian function. '^^ There- 
fore, on account of our variational results reported in 
Sec. nil Al which show that molecule exchanges are likely 
to occur at very low rate, it may be suggested that 
two-dimensional hydrogen presents some degree of an- 
harmonicity. 

In Fig. [8] (Top), we show DMC results similar to 
those enclosed in Fig. [7] but for larger distances and 
expressed in logarithmic scale in order to obtain the 
asymptote of gi . As one observes in there, the value 
of limr^oo gi in the unsymmetrized case tends obviously 
to zero, while for 4'jg i^ amounts to a small but finite 
value no = 6(1) • 10"'' . In the same figure (Bottom), 
we compare gi(r) at several densities and only for the 
symmetric wave function; we obtain no = 2(1) • 10""' 
and 8(1) • 10"^ at the density 0.067 and 0.056 A'^ , re- 
spectively. Apart from the fact that we obtain a small 
but finite condensate fraction for solid P-II2 in all the 
studied cases, we note that the value of no raises very 
abruptly in moving from equilibrium to densities close 
to the spinodal point (where P < 0). In a very recent 
work, we have analyzed the superfluid nature of solid ^He 
at zero temperature by means of ipjQ ^ In that work, 
we showed that the superfluid fraction of bulk solid ^He 
lies below 1 • 10~^, whereas a clear superfluid signal of 
Ps/p — 3.2(1) • 10"'^ appears in the presence of 1 % of 
vacancies. In the case of vacancies, we found that the 
condensate fraction increased by roughly a factor two 
with respect to that of the perfect crystal configuration. 
According to this outcome, a significant increase of no in 
our simulations might be identified to the appearance of 
superfluidity in the system. 

In oder to ensure that the signiflcant raise of no ob- 
served in our simulations is not due to partial melt- 
ing of H2, we have calculated the corresponding radial 
pair distribution function at the density 0.056 A~^ . In 
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FIG. 9: Top: Radial pair distribution function of two- 
dimensional H2 and He near the spinodal and freezing den- 
sities, respectively. Bottom: Structure factor S{k) of 2D H2 
at equilibrium density po and p — 0.056 A~^ . 



Fig. [9] (Top) , we report g{r) for hydrogen and compare it 
with the one obtained for two-dimensional ^He above its 
freezing point. Clearly, a typical solid pattern emerges for 
H2 . Moreover, in the same figure (Bottom) we also plot 
the structure factor S(k) for molecular hydrogen at the 
equilibrium density po = 0.067 A~^ and p = 0.056 A^^ ; 
in both cases marked peaks emerge at the reciprocal lat- 
tice vectors. Therefore, the significant variation of uq 
that we observe when the density decreases is not caused 
by possible structural instabilities affecting the solid. 



VI. DISCUSSION AND CONCLUSIONS 



ratio and quoted so isotopic quantum effects in hydrogen. 
Our results show that no stable liquid phase exists and 
therefore reducing one dimension with respect to bulk 
it is not enough to get the so-longly searched superfluid 
phase of H2 . 

Interestingly, Wiechert et al. have reported very re- 
cently on an experiment on molecular ortho-deuterium 
coadsorbed on graphite preplated by a layer of Kr, up 
to temperatures of ^ 1.5 Kj^ The authors of this work 
claim evidence for the existence of a reentrant D2 liq- 
uid at very low temperatures, based on their heat ca- 
pacity and neutron diffraction measurements. The sys- 
tem explored by Wiechert et al. can be fairly modeled 
by a monolayer. However, on account of our results for 
H2, the possibility of pure two-dimensional liquid deu- 
terium at zero temperature must be ruled out. On the 
expectance of new and more explanatory experiments, we 
may point that, assuming that thermal effects are prac- 
tically negligible, the role of the interactions between the 
deuterium molecules and the atoms of the substrate are 
the ones of relevance. Certainly, TurnbuU and Boninsegni 
have already addressed recent work on this direction by 
means of the Path Integral Monte Carlo (PIMC) method 
and simple interaction modelsi^'^ Further improvement 
on the modeling of coadsorbed systems, putting especial 
emphasis on the description of the interactions and the 
effect of corrugation with the substrate, may open new 
and challenging venues for the realization of superfluidity 
in P-II2 and 0-D2 systems.— 

At the variational level, we have analyzed the quality 
of three different symmetrized trial wave functions based 
on the Nosanow-Jastrow model in describing 2D solid 
molecular hydrogen. We have shown that the recently 
proposed symmetrized wave function used to describe the 
supersolid also characterizes hydrogen satisfactorily. By 
using that wave function, we have studied the behavior 
of the one-body density matrix of solid P-H2 with density 
and predicted that the system could become superfluid 
at very dilute densities (where P < 0). Further work is 
being carried out to estimate the superfluid density in the 
negative pressure region trying to confirm the signature 
observed in the one-body density matrix. 



To summarize, in this work we have studied two- 
dimensional P-II2 and 0-D2 at zero-temperature and low 
pressures, with the diffusion Monte Carlo method and 
the Silvera-Goldman semi-empirical pair interaction. We 
have assessed several energetic and structural properties 
of both systems, like the total and kinetic energy per par- 
ticle, radial pair distribution function, and Lindemann's 
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